1 확진자 데이터

COVID 팩키지에서 한국에서 코로나19 검사결과 확진된 시계열 데이터를 얻는다. 시계열 데이터 예측을 위한 기초 데이터를 준비한다.

library(tidyverse)
library(COVID19)
library(timetk)

korea_dat <- covid19("KOR", level = 1)
We have invested a lot of time and effort in creating COVID-19 Data Hub, please cite the following when using it:

  Guidotti, E., Ardia, D., (2020), "COVID-19 Data Hub", Journal of Open
  Source Software 5(51):2376, doi: 10.21105/joss.02376.

A BibTeX entry for LaTeX users is

  @Article{,
    title = {COVID-19 Data Hub},
    year = {2020},
    doi = {10.21105/joss.02376},
    author = {Emanuele Guidotti and David Ardia},
    journal = {Journal of Open Source Software},
    volume = {5},
    number = {51},
    pages = {2376},
  }

To retrieve citation and metadata of the data sources see ?covid19cite. To hide this message use 'verbose = FALSE'.
korea_df <- korea_dat %>% 
  ungroup() %>% 
  dplyr::select(날짜=date, 누적검사자=tests, 누적확진자=confirmed, 누적회복자=recovered, 누적사망자=deaths) %>% 
  timetk::pad_by_time(.date_var = 날짜, .by = "day", .pad_value = NA) %>% 
  mutate(누적검사자 = timetk::ts_impute_vec(누적검사자, period = 1) %>%  ceiling(.),
         누적회복자 = ifelse(is.na(누적회복자), 0, 누적회복자),
         누적사망자 = ifelse(is.na(누적사망자), 0, 누적사망자))

korea_daily_df <- korea_df %>% 
  mutate(검사자 = 누적검사자 - lag(누적검사자, n = 1L),
         확진자 = 누적확진자 - lag(누적확진자, n = 1L),  
         회복자 = 누적회복자 - lag(누적회복자, n = 1L),  
         사망자 = 누적사망자 - lag(누적사망자, n = 1L) ) %>% 
  drop_na() %>% 
  select(-contains("누적"))

confirmed_daily_tbl <- korea_daily_df %>% 
  select(날짜, 확진자) %>% 
  timetk::pad_by_time(.date_var = 날짜, .by="day", .pad_value = NA)

confirmed_daily_tbl 
# A tibble: 339 x 2
   날짜       확진자
   <date>      <int>
 1 2020-01-23      0
 2 2020-01-24      1
 3 2020-01-25      0
 4 2020-01-26      1
 5 2020-01-27      1
 6 2020-01-28      0
 7 2020-01-29      0
 8 2020-01-30      0
 9 2020-01-31      7
10 2020-02-01      1
# ... with 329 more rows

2 시각화

timetk 팩키지 plot_time_series() 함수를 사용해서 확진자 추세를 시각화한다.

confirmed_daily_tbl  %>% 
  plot_time_series(.date_var   = 날짜, 
                   .value      = 확진자, 
                   .line_color = "#ff0000",
                   .title      = "코로나19 일별 확진자 추세")

3 모형 데이터

3.1 전체 데이터

confirmed_daily_tbl 데이터를 history_tbl 데이터프레임으로 저장시켜 명확히 의미를 부여한다. timetk 팩키지 future_frame() 함수로 예측할 기간(“1달”)으로 데이터프레임을 확장하여 관측된 시계열과 예측 시계열을 결합시킨 전체 데이터프레임 데이터를 준비한다. forecast_tbl 로 예측할 데이터프레임을 따로 떼어둔다.

# 전체 시계열 데이터 ----
full_tbl <- confirmed_daily_tbl %>% 
  timetk::future_frame(.length_out = "1 month",
                       .bind_data  = TRUE)
# 관측 시계열 데이터 ----
history_tbl <- full_tbl %>% 
  filter(!is.na(확진자))

# 예측 시계열 데이터 ----
forecast_tbl <- full_tbl %>% 
  filter(is.na(확진자))

full_tbl %>% 
  plot_time_series(.date_var = 날짜, 
                   .value    = 확진자,
                   .title    = "코로나19 일별 확진자수")

3.2 훈련/시험 데이터

time_series_split() 함수를 사용해서 훈련/시험 표본으로 데이터를 구분한다. assess = 28을 지정하여 향후 28일 후를 예측하는 모형을 개발한다.

splits <- history_tbl %>%
    time_series_split(날짜, assess = 28, cumulative = TRUE)
    # time_series_split(날짜, initial = "3 month", assess = 28, cumulative = FALSE)

splits
<Analysis/Assess/Total>
<311/28/339>

상기 훈련/시험 데이터 구분한 후 시각화를 통해 확인해보자. timetk::plot_time_series_cv_plan() 함수를 사용한다.

splits %>%
    tk_time_series_cv_plan() %>%
    timetk::plot_time_series_cv_plan(날짜, 확진자)

4 피처 공학

tidymodels 생태계의 recipe 팩키지 recipe() 함수를 사용해서 피쳐 공학(feature engineering) 작업을 수행한다. 추후, 다양한 피처를 작업할 예정이라 대략적인 틀만 잡아둔다.

library(tidymodels)

recipe_spec <- recipes::recipe(확진자 ~ ., data = training(splits))

recipe_spec %>% prep() %>%  juice()
# A tibble: 311 x 2
   날짜       확진자
   <date>      <int>
 1 2020-01-23      0
 2 2020-01-24      1
 3 2020-01-25      0
 4 2020-01-26      1
 5 2020-01-27      1
 6 2020-01-28      0
 7 2020-01-29      0
 8 2020-01-30      0
 9 2020-01-31      7
10 2020-02-01      1
# ... with 301 more rows

5 모형적합

parsnip 팩키지의 linear_reg() 를 통해 선형회귀(?) 모형을 시계열 데이터에 날짜를 독립변수로 넣어 말도 되지 않지만 가장 단순하게 예측 모형을 적합시킨다. 이를 위해서 workflow() 를 사용하는데 기본적으로 피처 공학의 recipe_spec이 필요하고 model_spec의 모형이 필요하다.

그리고 fit() 함수를 사용해서 적합을 시킨다.

model_spec <- linear_reg(
    mode = "regression") %>% 
  set_engine("lm")

wkfl_fit_lm <- workflow() %>% 
  add_recipe(recipe_spec) %>% 
  add_model(model_spec) %>% 
  fit(training(splits))

wkfl_fit_lm
== Workflow [trained] ==========================================================
Preprocessor: Recipe
Model: linear_reg()

-- Preprocessor ----------------------------------------------------------------
0 Recipe Steps

-- Model -----------------------------------------------------------------------

Call:
stats::lm(formula = ..y ~ ., data = data)

Coefficients:
(Intercept)         날짜  
 -5725.9014       0.3164  

6 모형 평가

workflow 객체를 modeltime 팩키지 modeltime_table() 함수에 넣어 객체로 만든 후에 modeltime_accuracy() 함수를 사용해서 모형 성능을 파악한다.

library(modeltime)

model_tbl <- modeltime_table(
    wkfl_fit_lm
)

model_tbl %>% 
  modeltime::modeltime_accuracy(testing(splits))
# A tibble: 1 x 9
  .model_id .model_desc .type   mae  mape  mase smape  rmse   rsq
      <int> <chr>       <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1         1 LM          Test   661.  78.4  6.48  130.  701. 0.802

6.1 시각화

시험데이터를 통해 모형으로 예측한 값을 시각화한다. \(R^2\) 값도 80%가 나오는데 시각화를 하니 시험데이터에 대한 적합이 데이터와 동떨어진 것을 확인할 수 있다.

calibration_tbl <- model_tbl %>% 
  modeltime_calibrate(
    new_data = testing(splits)
  )

calibration_tbl %>% 
    modeltime_forecast(
        new_data      = testing(splits),
        actual_data   = history_tbl,
        conf_interval = 0.10
    ) %>%
    plot_modeltime_forecast(
        .legend_max_width = 60,
        .legend_show = FALSE,
        .conf_interval_show = TRUE,
        .conf_interval_alpha = 0.20,
        .conf_interval_fill = "lightblue",
        .title = "코로나19 확진자 1개월 예측 - 선형회귀"
    )

7 확진자 예측

앞선 모형은 history_tbl 을 훈련/시험 데이터에 대해 적합을 시킨 것이라 … 이제 시간을 확대하여 modeltime_refit() 함수를 사용해서 모형을 전체 데이터에 대해 다시 적합시킨다.

refit_tbl <- calibration_tbl %>% 
  modeltime_refit(data = history_tbl)

마지막으로 앞서 구축된 모형을 바탕으로 현재까지 입수된 데이터를 바탕으로 한달 후를 예측한다.

refit_tbl %>% 
  modeltime_forecast(
    new_data    = forecast_tbl,
    actual_data = history_tbl,
    conf_interval = 0.5
  )  %>%
    plot_modeltime_forecast(
        .legend_max_width = 25,
        .conf_interval_fill = "lightblue",
        .conf_interval_alpha = 0.7,
        .interactive = FALSE
    )

8 모형과 데이터 저장

마지막 단계로 데이터와 모형을 저장하여 다음 단계로 나가기 위한 작업을 수행한다.

# 데이터 ----
full_tbl %>% 
  write_rds("data/full_tbl.rds")

# 모형 ----
wkfl_fit_lm %>% 
  write_rds("data/wkfl_fit_lm.rds")
 

데이터 과학자 이광춘 저작

kwangchun.lee.7@gmail.com